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Abstract 

We report our rederivation of the equations of motion for relativistic compact binaries through 
the third-and-a-half post-Newtonian (3.5 PN) order approximation to general relativity using the 
strong field point particle limit to describe self-gravitating stars instead of the Dirac delta func- 
tional. The computation is done in harmonic coordinates. Our equations of motion describe the 
orbital motion of the binary consisting of spherically symmetric non-rotating stars. The resulting 
equations of motion fully agree with the 3.5 PN equations of motion derived in the previous works. 
We also show that the locally defined energy of the star has a simple relation with its mass up to 
the 3.5 PN order. 
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I. INTRODUCTION 



A relativistic compact binary (e.g., a neutron star binary) loses its orbital angular mo- 
mentum by emitting gravitational waves and coalesces in the end. Such a system is a 
promising source for the gravitational wave detectors such as CLIO, GEO600, LCGT, LIGO, 
TAMA300, and VIRGO jj. However, even with those advanced detectors, a direct detec- 
tion is not easy. Indeed, because the amplitude of gravitational wave from such a source is 
expected to be tiny at the Earth compared to the detectors' noises, an efficient detection 
method is required. One method widely used is matched filtering. When using this tech- 
nique, it is known that the more accurately we know the shape of the signal, the larger the 
signal to noise ratio becomes. This in turn means that it is desirable to know the dynamics of 
the binary accurately when one hopes to increase the number of the detectable events. This 
purpose can be achieved by using higher order post-Newtonian (PN) equations of motion 
for two point particles, because up to the last several orbits before the coalescence the bi- 
nary components have moderately slow orbital velocities and are affected negligibly by tidal 
effects . In fact, it is suggested that at least a third order post-Newtonian correction 
in the equations of motion are necessary for extraction of astronomical information of the 
sources |4|. To obtain higher order post-Newtonian corrections, mainly three methods are 
employed in the literature. 

In the first method, the Arnowitt-Deser-Misner (ADM) Hamiltonian is derived by means 
of a direct PN iteration of the Einstein equations 
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lia|. Indeed, 

the ADM Hamiltonian in the ADM Transverse- Tarceless (TT) gauge is completed up to the 



3.5 PN order inclusively 
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In the second method, one assumes that the energy and 
the angular momentum fluxes at infinity are balanced by the corresponding loses of those 
in the binary orbital motion. Here the known PN expressions for the energy and angular 



momentum fluxes at infinity are used. With this second method, the n PN (n = 2.5, 3.5, 
and 4.5) order corrections to the Newtonian equations of motion are derived [l6, 17, fla. fioj] - 



Finally in the third met 
the Einstein equations 



rodjPN equations of motion are derived by a direct 
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istein equations [20|, l2ll. |22|. |23 
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iteration of 



35|, 136|, 137|, 



45j for reviews). Wit 



of motion are completed up to the 3.5 PN order inclusively 



i this method the equations 
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38]. 



This paper reports a rederivation of the 3.5 PN correction to the Newtonian acceleration 
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for a spherically symmetric non-rotating self-gravitating star using the same method as that 



in our previous papers 30 



the 3 PN order (see 



33 



35 



36 



371 ] which derived the corrections up to and including 



45j| for a review). We use the third method mentioned in the previous 



paragraph. However our method is different in the several aspects from those of the previous 



works 



34. 



381 ] which also used the third method to derive the 3.5 PN order corrections. Here 
we briefly mention the two most important points. 

The first point is regarding the way to describe the binary components as point parti- 
cles. Among the previous works, Damour, Blanchet, Schafer and their collaborators have 
used Dirac delta functionals to describe binary stars as point particles. However, divergent 
integrals appear when using Dirac delta functionals in general relativity. Then to regularize 



those divergences, for example, the work 



successfully used the Hadamard Partie Finie 



(HPF) regularization to derive the 2.5 PN equations of motion in harmonic coordinates. 
Blanchet and Faye have developed the generalized HPF regularization 46j in a Lorentz 
invariant manner 47j and derived the 3 PN correction in the same gauge except for one 
and only one numerical coefficient (denoted by A) which could not be determined within 



their method 



29 



321 ] . Interestingly, two works [4i 



have shown that this coefficient 



A actually corresponds to one of the undetermined coefficients (u; sta tic) reported previously 
by Jaranowski and Schafer in their derivation of the 3 PN ADM Hamiltonian in the ADM 
Transverse- Traceless (TT) gauge jll, l^Q]- Indeed these latter works also used Dirac delta 
functionals and the HPF. It has become clear that the HPF is not an appropriate method 
for regularization of the divergent integrals at the 3 PN order. Damour, Jaranowski and 
Schafer Jjj then used the dimensional regularization to derive the 3 PN ADM Hamiltonian 
in the ADMTT gauge and they finally completed the 3 PN correction that contains no un- 
determined coefficient. The 3 PN equations of motion in harmonic coordinates were later 
derived using the combination of the HPF and the dimensional regularization 50(. Their 
result physically agree with the result of pjj]. Incidentally, there found no trouble when 



using the HPF at the 3.5 PN order and the works 10 
order using the HPF. See J] for a review. 



38] derived the corrections at that 



271 ] proposed the 



Now what is our method to achieve a point particle limit? Futamase 
strong field point particle limit. In this limit the binary star is first described as an extended 
object and then its radius is taken to zero in a specific manner which will be explained in 
the Section II below. By this limit we obtain a point particle with strong internal gravity 
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that is appropriate as a model of a relativistic compact star. With this limit and the method 



mentioned below we have successfully derived the 3 PN equations o 
no undetermined coefficient 



35 



361 ] and confirmed the earlier result 



motion that contains 



The second aspect of our method that is different from the others is the way to derive 
equations of motion. In our method, using the local energy momentum conservation law, we 
compute the force acting on the star by computing the gravitational energy momentum flux 
going through a suitably defined surface around the star. This idea was used by Einstein, 



Infeld and Hoffmann in their derivation of the 1 PN equations of motion 



other hand, it was assumed in some works ( 28( at the 2.5 PN order, 



order, and 



3. 



29 
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2l|. On the 



32J at the 3 PN 



at the 3.5 PN order) that the binary star follows a "geodesic" of spacetime. 



However for an equal mass binary, the metric components diver ge a t the position of the star 



29|, 



32 



the divergence 



when it is described as a point particle. In those works 
was regularized by the HPF (at the 2.5 and 3.5 PN order), or the combination of the HPF 
and the dimensional regularization (at the 3 PN order). Our results up to the 3 PN order 
inclusively have fully agreed with these results. This agreement has confirmed that the star 
follows the regularized geodesies at least through the 3 PN order. One of the aims of the 
current paper is to extend it to the 3.5 PN order. (In this respect, Fukumoto et al 51 ] 
showed that the star should follows the geodesic of the "smooth part of the full metric" 
where the smoothing is done by means of the surface integral approach. However, this work 
showed it for the case of an extreme mass ratio binary.) 

This paper is organized as follows. The next section explains our method of deriving the 
equations of motion, sketching each steps in our method with emphases on the strong field 
point particle limit, the surface integral approach and the scalings of the stress energy tensor 
of the matter assumed in our approach. Sec. II I II explains how to derive the PN gravitational 
field and the PN equations of motion. After showing the formal structure of the 3.5 PN 
equations of motion in Sec. IIVI the 3.5 PN gravitational field in harmonic coordinates is 
derived in Sec. |V] Since we shall not use a Dirac delta functional nor assume any specific 
form of a momentum velocity relation, we have to derive the relation between the star's 
mass and its energy and the one between the star's three momentum and its three velocity. 
Those relations at the 3.5 PN order are reported in the sections IVII and [VIII We show the 
resulting 3.5 PN equations of motion in Sec. IVIIII 

Throughout this paper, we will use harmonic coordinates and the unit of c = 1 = G 
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where c is the velocity of light and G is the Newton's gravitational constant. A tensor with 
alphabetical indexes such as x\ denotes a Euclidean three- vector. We raise or lower its 
indexes with a Kronecker delta. For an object with Greek indexes such as x M , its indexes 
are raised or lowered with a flat Minkowskian metric. We shall call our previous papers 



3fJ, |33j, |36( Paper I, II, and III, respectively. 



II. KEY IDEAS IN OUR DERIVATION 



In this section we briefly review our method to obtain post-Newtonian equations of mo- 
tion, listing the procedures in our method with emphases on the strong field point particle 
limit, the surface integral approach and the associated scalings of the stress energy tensor 



of the matter. See 



27 



30 



33 



42 



44 



|52J,[53J,[54J for more details. 
We first introduce an adimensional parameter e which represents the slowness of the star's 
orbital velocity v z orh . 



V ° Th ~ dt ~ € dr 



where z l is the star's representative point and we assume v % = dx 1 1 dr of order unity. The 



27|, 



44j . The post-Newtonian 



time coordinate r is called the Newtonian dynamical time 
scaling implies that the square modulus of the orbital velocity is of order of the inter-body 
gravity fh/r where fh and r are the star's mass and the orbital radius, respectively. In 
terms of e, this means that the mass of the star fh is of order of e 2 3fJ, |33j, |36j. In what 
follows, m (without a tilde) denotes the star's mass that is of order unity. The smallness 
of the orbital velocity and the inter-body gravity are represented in terms of e. Thus, e is 
a post-Newtonian expansion parameter. Henceforth we call the (r, x l ) coordinate the near 
zone coordinate. 

The typical magnitude of the velocity of the system far away from the binary is the 
velocity of the gravitational wave emitted by the binary, or the velocity of light c. Hence 



there must be different coordinates appropriate in the far zone [55|, [56]. The near zone 
coordinate is useful within, say, one wavelength of the gravitational waves (of order ~ cx 
(orbital period of the binary)) from the binary system center of the mass. 

Now, the next subsection sketches the principal steps of our derivation of post-Newtonian 
equations of motion. 
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A. Principal steps of our derivation 



Suppose we have the n PN metric and the n PN equations of motion. To derive the n + 1 
PN equations of motion, we first solve the Einstein equations to the n + 1 PN order and 
then compute the gravitational force acting on the binary star. When solving the Einstein 
equations, we use the following techniques. 

• Gauge: We use harmonic coordinates to solve the Einstein equations. 

• Einstein equations: We solve the harmonically relaxed Einstein equations where we 
use the Minkowskian wave operator for computational simplicity and where we rewrite 



the Einstein equations as a set of retarded integrals in the same manner as in 



Boundary condition at infinity: We take the no-incoming radiation condition at the 



Minkowskian past null-infinity 



41 



Other previous works 34j, |38( take essentially the same approaches on all of the above points. 



Now it is well-known that the post-Newtonian a ppr oximation breaks down far away from 



the binary at high post-Newtonian orders (e.g., [59|]). To obtain the metric, it is thus 
convenient to divide the spacetime into the near zone and the far zone. The near zone is a 
time-like world tube surrounding the binary stars. The center of the near zone is defined (if 
necessary) as the center of the mass of the binary. The intersection between the r =constant 
spatial hypersurface and the world tube becomes a spatial 3-sphere. The radius of the sphere 
is of order of the largest wavelength of the gravitational wave emitted by the binary due 
to its orbital motion. Its mathematical definition shall be given later in Sec. IIII Al Then 
the far zone is the outside of the near zone. The domain of the integration of the retarded 
integrals is now divided into that in the near zone and that in the far zone. We compute 
two contributions separately. 

In the near zone we can use the post-Newtonian approximation safely. We use the 
following method to obtain the near zone contribution. 

• Near zone contribution: To obtain the near zone metric, we first expand the retarded 
integrals in e and change the domain of the integration from the Minkowskian null-cone 
to the r = constant spatial hypersurface. This hypersurface is denoted as NZ. 
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• Body zone: The region NZ is divided into three regions. Two of the three are spheres 
called the body zones each of which surrounds each of the binary star. The remaining 
is the region outside of the two body zones. 

• Multipole expansion: In each body zone, the integrals are evaluated by a multipole 
expansion. 

In the multipole expansion, the multipole moments are defined as the volume integrals on 
the r = constant spatial hypersurface whose integrands include the matter's stress energy 
tensor plus the gravitational stress energy tensor. In the usual post-Newtonian approxima- 
tion, one assumes that the gravitational field is everywhere weak. Since the gravitational 
wave detectors search for relativistic compact binaries that have strong internal gravity, we 
would like to have a method where we use the post-Newtonian approximation only outside of 
;he stars. One procedure which we adopt here to solve this problem is proposed by Futamase 



27] . 



• Strong field point particle limit [27J]: With this limit the star's internal gravity can 
be assumed to be strong. Namely if the companion star were absent, the star's mass 
would become the ADM mass. We shall explain this limit in more detail in the next 
subsection. 

As in the Newtonian dynamics, the multipole moments are defined with respect to some 



reference point 



60]. 



• Star's representative point: The representative point of the star (e.g., the center of 
the mass of the star) is defined by setting the star's dipole moment appropriately in 
the same manner as in the Newtonian dynamics. The star's multipole moments are 
defined with respect to this representative point. 

As is mentioned in the introduction, this paper studies an intrinsically spherically sym- 
metric star. As the multipole moments above are defined as volume integrals on the r = 
constant spatial hypersurface, we cannot obtain a spherically symmetric star by simply mak- 
ing the multipole moments vanish. This is because the Lorentz contraction makes such a 
star acquire apparent multipole moments when an observer sees the star moving with respect 
to the observer. We define the star's intrinsic multipole moments in the generalized Fermi 
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normal coordinate 61] and set those zero to obtain an intrinsically spherically symmetric 
star. 

Now, in the far zone, it is well-known that a simple post-Newtonian expansion gives 



divergent integrals (e.g., 59j). To solve this problem, we use the following method. 



Far zone contribution: we use the Direct Integration of the Relaxed Einstein equations 
(DIRE) method to evaluate the far zone contribution. 



This method was proposed by Will and his collaborators 3l|, |34|, [62j. This method assumes 



that the binary is sufficiently stationary in the long past. On the other hand, the work [38] 
used the multipolar post-Minkowskian (MPM) approach [63] to obtain the far zone metric. 
The asymptotic matching is done between the PN inner metric and the MPM outer metric 



64]. 



Having the metric components to the n + 1 PN order, we now derive the n + 1 PN 
equations of motion for relativistic compact binaries. The following steps are used only in 
our method. 

• Surface integral approach: Using the local conservation law of the stress energy, we 
obtain expressions for the time derivative of the four momentum in terms of the surface 
integrals over the star's body zone surface. The integrands shall be the Landau-Lifshitz 
pseudo-tensor. Thus, the surface integrals give the net flux going through the star's 
body zone surface and amount to the gravitational force acting on the star. 

• Four momentum - mass and velocity relation: As is written in the previous point, we 
compute the time derivative of the four momentum. To obtain equations of motion, 
we need to compute the time derivative of the coordinate velocity. In our paper the 
star's four momentum is defined as the volume integral of the sum of the matter's 
and gravitational stress energy tensors. It is non-trivial that the four momentum is 



proportional to the four velocity 60j. We develop a method to derive the expression 



for the four momentum in terms of the stars' masses and their coordinate velocities. 

Gravitational energy momentum flux: Using the metric components (the solution of 
the Einstein equations to the n + 1 PN order), we evaluate the Landau-Lifshitz pseudo- 
tensor to the n + 1 PN order. 



S 



Surface integrals: We evaluate the surface integrals to the n + 1 PN order and obtain 
the n+1 PN equations of motion. 



The previous works 



31 



34[ use the volume integral approach where they assume a geodesic, 



multiply it by the conserved baryon density, and integrate it over the star on the time 



constant hypersurface. Nissanke and Blanchet 38] also assume a geodesic, but in their 
method they substitute the metric components into the geodesic and regularize divergence 
due to their use of Dirac delta functionals. 

In the following subsections, we shall explain in more detail the strong field point particle 
limit, the surface integral approach, and the associated scalings of the stress energy tensor 
components. Other features listed above shall be explained in Sec. IIHI We apply the 
procedures listed in this section to the derivation of the 3.5 PN order equations of motion 
starting from Sec. HVl 



B. Strong Field Point Particle Limit 

One would describe a star as a point particle by making the radius of the star zero. 
However, by this procedure, one would obtain a black hole rather than a point particle. 
The radius of the star cannot be made smaller than of order of its mass (in the unit of 
G = 1 = c). Futamase 



271 ] proposed that we may obtain a point particle model for the star 
by taking a limit where both the radius of the star and its mass shrink at the same rate. 
Since the post-Newtonian scaling implies that the star's mass m is 0(e 2 ), we assume that 
its radius is also 0(e 2 ). The point particle limit is now achieved in the limit where e goes to 
zero. 

The usual post-Newtonian approximation assumes the gravitational field is everywhere 
weak even inside the stars. On the other hand, with Futamase's procedures we obtain a 
point particle with finite internal gravity, since a typical magnitude of the self-gravitational 
field (the mass over the radius) of the point particle achieved by this procedure is finite 
irrespective of e. For this reason, this limit is called the strong field point particle limit. 
Note that the inter-body gravity (the mass over the orbital separation) is 0(e 2 ) and the 
post-Newtonian approximation applies for the aorbital motion. With this limit, we can 
derive post-Newtonian equations of motion for a binary star whose component stars have 
strong internal gravity. 



C. Surface Integral Approach and Body Zone 

One way to derive equations of motion of a star is a surface integral approach where 
we compute the total gravitational energy momentum flux going through a surface around 
the star. When we let the radius of the surface shrink to zero and when the radius of the 
star goes to zero faster than that of the surface in the point particle limit, we would then 
obtain equations of motion for the star as a point particle. For this purpose, we introduce 
two body zones B A for the stars A = 1,2 as B A = — za{t)\ < €Ra}- Here z 1 a (t) is a 

representative point of the star A, e.g., the center of the mass of the star A. We shall later 
fix z A by specifying the star's dipole moment as in the Newtonian dynamics. The body zone 
radius eR A are much smaller than the orbital separation but is larger than the radius of the 
star for any e (Recall that the radius of the star decreases proportionally to e 2 in the strong 
field point particle limit while the body zone radius does so proportionally to e). Note that 
the two body zones do not overlap each other. Finally, Ra are constant, i.e., (IRa/cIt = 0. 
Other than these conditions, Ra are arbitrary. 

Now the above scalings of the masses and the radii of the stars motivate us to introduce 
a body zone coordinate for the star A as (r, or A ) where a~ A = e~ 2 (x % — z A (r)). The scalings 
of the body zones and the body zone coordinates give us a situation where in the body zone 
coordinate A the star A does not shrink as e — > while the boundary of the body zone 
expands to infinity. Thus, it is appropriate to define the star's characteristic quantities such 
as its mass using the body zone coordinate. Moreover, since the body zone boundary dB A 
is far away from the surface of the star A (in its body zone coordinate), we can evaluate 
explicitly the gravitational energy momentum flux on dB A using the post-Newtonian gravi- 
tational field. After evaluating the surface integrals, we make the body zone shrink to derive 
the equations of motion for the compact star. 

Possible effects of the internal structures of the compact stars are coded in the multipole 
moments of the stars. These moments in turn appear in the gravitational energy momentum 
flux and would affect the orbital motion. However, in this paper we shall concentrate on 
spherically symmetric stars and ignore those multipole moments. 



10 



D. Scalings of the Matter Stress Energy Tensor 



The scalings of the radii and the masses of the stars indicate that the matter density 
is 0(e~ 4 ) in the (t, x l ) coordinate (or e~ 2 in the near zone coordinate (r, x 1 )). We further 
assume that the internal time scale of the star is comparable to that of the binary orbital 
motion and is 0(e). Namely, in this paper we assume that the star is pressure supported 
and non-rotating, although an extension to rapidly rotating stars is straightforward 27]. 

In terms of the stress energy tensor of the matter T^ u (or the source terms of the harmon- 
ically relaxed Einstein equations in Eq. (T5]) below), these scalings imply T TT = 0(e -2 ), 
T T ~ = 0(e~ 4 ), T^- = 0(e~ 8 ) where the underlined indexes mean that for any tensor S\ 



271 ] . Because a 



S l = e~ 2 S % and reminds the scaling of the body zone spatial coordinate 
star moves in a gravitational field, the tensor components in the body zone coordinate are 
different from those in the near zone coordinate. Let us denote the matter's stress energy 
tensor in the body zone coordinate (r, a A ) by T A and that in the near zone coordinate 
(r, x l = z' l A + e~ 2 a A ) by T^ z . Transforming to the near zone coordinate we obtain 271 ] 



rpTT rprr 

1 NZ — 1 A i 



T N l z = eT/ + v\T\\ 

T^ z = e*T A t + 2e 2 v A rr A )T + v A v A T A - 



Hence the stress energy tensor components of the matter in the near zone coordinate varies 



with respect to e in the same way as those in the body zone coordinate, or in short, T 
0(e- 2 ), Tjf z = 0( e - 4 ), T% = 0(e~ 8 ) £]. 



TT 

NZ 



III. MATHEMATICAL FORMULATION 

Based on the idea explained in the previous section, this section formulates our method 
to derive the equations of motion. We first explain how to solve the Einstein equations, then 
show how we achieve the surface integral approach. 

A. Field Equation 

In the surface integral approach, we need to compute the gravitational field near the body 
zone boundary where the field is well described by the post-Newtonian approximation and 
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slightly deviates from the flat metric. We define a deviation field h^ u as 

hT = rf v - V^gg' u , (i) 

where if" = diag(— e 2 , 1, 1, 1) is the flat metric in the near zone coordinate (r, x l ) and g is 
the determinant of the metric. The indexes of h^ u are raised or lowered by the flat metric. 

Now we impose a harmonic coordinate condition = where the comma denotes 

a partial derivative. In the harmonic gauge, we can recast the Einstein equations into a 
relaxed form, 

UN™ = -16ttA^, (2) 

where □ = rj^d^du is the flat spacetime d'Alembertian. The source term A Miy of the relaxed 
Einstein equations consists of two pseudo-tensors. The first is the sum of the stress energy 



tensor of the stars denoted by T^ v and the Landau-Lifshitz pseudo-tensor t^ L 65j. The 
second arises due to our use of the flat spacetime d'Alembertian instead of the curved 
spacetime one. The explicit expressions are 

A^ = 6^ + x^,*?, (3) 

&» = {-g)(T' a ' + t% i ), (4) 
x ^af3 = _( h av h f3n _ ^h^). (5) 

The harmonic condition on h^ v implies a local energy momentum conservation law; 

h»\ v = 0. (6) 

Note that x^ ual3 ,ap itself is divergence free, i.e., X^ a ^,a/3u = 0. 

We can formally rewrite the relaxed Einstein equations as retarded integrals; 



h^(r, = 4 / cfy " V ,J J"" (7) 

Jc(r,x k ) \ x ~ y\ 

where C(r,x k ) means the past light cone emanating from the event (r, x k ). We have as- 
sumed no homogeneous solution of the relaxed Einstein equations. It is well-known that 
this condition can be deduced from the Minkowskian no-incoming radiation condition (See, 



e.g., the section 92 of 58] or the section 6 of [41]). 



We solve the Einstein equations as follows. First we split the domain of the integration 
into two zones: the near zone and the far zone. The near zone is the neighborhood of 
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the gravitational wave source where the wave character of the gravitational radiation is not 
manifest. In this paper, as in our previous papers, we define the near zone as a time-like 
world tube surrounding the binary stars. The center of the near zone is defined (if necessary) 
as the center of the mass of the binary. The intersection between the r =constant spatial 
hypersurface and the world tube becomes a spatial 3-sphere. Mathematically, denoting the 
harmonic coordinate distance from the center by r, the near zone is defined as r < IZ/e 
where IZ/e is of order of the largest wavelength of the gravitational wave emitted by the 
binary due to its orbital motion. We assume TZ/e sufficiently large so that the near zone 
covers the binary stars. Finally, 1Z is constant in time, i.e., dlZ/dr = 0. Otherwise 1Z is 
arbitrary. The e -1 scaling of the near zone radius is derived from the e dependence of the 
wavelength of the gravitational wave emitted by the binary. The outside of the near zone is 
the far zone where the retardation effect of the field is crucial. 

For the retarded integrals in the far zone, we evaluate it using the Direct Integration of 
the Relaxed Einstein equations (DIRE) method. This method was proposed by Will and his 



collaborators 



31 



34 



621 ] . DIRE directly and nicely fits into our formalism since it utilizes 



the relaxed Einstein equations in the harmonic gauge. Using DIRE, Pati and Will 31 ] 
showed that the far zone contribution to the near zone field affects (physically) the orbital 
motion starting at the 4 PN order In fact, this result was first obtained by Blanchet and 
Damour [66j. Although we do not show our explicit computation in this paper, we have 
followed the DIRE method and checked that the far zone contribution does not affect the 
equation of motion through the 3.5 PN order. We shall thus focus our attention on the near 
zone contribution to the near zone field and neglect all the far zone contributions to the near 
zone field. 



B. Near Zone Contribution 



For the near zone contribution, we first make retardation expansion and change the 
domain of the integration to a r = constant spatial hypersurface 

h^rj) =4£^ (§X / d 3 y\x-yr^ z (r,y k ;e). (8) 

n=0 ' ^ ' JNZ 

where NZ denotes the near zone and we attach the subscript NZ to A M1/ to clarify that they 
are quantities in the near zone. Note that the above integral depends on the arbitrary length 
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1Z in general. The cancellation between the 1Z dependent terms in the far zone contribution 
and those in the near zone contribution through all the post-Newtonian order was shown 



by Pati and Will [31|- This paper uses their method, and hence we safely neglect all the 1Z 
dependent terms other than those where 7Z appears in arguments of logarithms. This is to 
make the arguments of all possible logarithmic terms adimensional. 

Second we split the integral into two parts; contribution from the body zone B = B1UB2, 
and from elsewhere, NZ/B. We thus evaluate the following two types of integrals 

^ = £ (hZ + h% z/Bn ) , (9) 



n=0 



'Vz/ fl „- 4 „, \ 9t ) J NZ/B dV \3-y\i-- {U> 
where ?a = x — Za- We shall deal with these two contributions successively in the followings. 



1. Body Zone Contribution 

As for the body zone contribution, we make multipole expansion using the scaling of the 
integrand (A^ z ) in the body zone. For example, the n = part in Eq. (ITU1) . h^ n=Q) gives 

/ P r n fc rr-k oT<kl>^k r l e T<klm>„k JL „m 

ttt .4 / A 1 ^2 -^A'A 1 4^fA ' A' A 1 ,6° j A 'a'a'A 

n Bn=o - 4e I 77 + e _ 7 h e 77 h e 



A=l,2 
.12 



+0(e 12 ), (12) 

(pi TfcvA; o j<kl>i k „/ \ 

^ + ^+^ \/ AA )+0^°), (13) 

_ ^ r A lr A J 

( yij y ki 'j r k q y<kl>ij h r l rr 7 <klm>ij h J, r m 

hV — A 2 I ZA _L ^ A A _l ^4 JZ; A r A r A , ^.6 0Z; A r A r A r A 

^Bn=0 " 4e 2^ 1 rA +£ r 3 + e 2r 5 + e 2r 7 



A=l,2 
.10 



+0(e lu ), (14) 

where = |r^|. The quantity with <> denotes a symmetric and tracefree (STF) operation 
on the indexes between the brackets. To derive the 3.5 PN equations of motion, we need 
h TT up to 0(e n ) and W* up to 0(e 9 ). 
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In the above equations we denned the multipole moments of the star A as 

Jf =e 2 / d 3 a A A% z af, (15) 



JB A 

J A ll = e 4 d 3 a A A T } ^ z af L , ( 16 ) 
JB A 

Z? ij = e 8 [ d s a A A% z o%, (17) 
Jb a 

where the capital index denotes a set of collective indexes, /; = and = a A L ar A 2 ---a A . 

Then P\ = J^ , Z)^ = I A , P^ 1 = J^ 1 . We simply call the four momentum of the star A, 
P\ the (three) momentum, and P\ the energy. Also we call D A the dipole moment and I A 
the quadrupole moment. 

Then we transform these moments into more convenient forms using the conservation law 
Eq. ([6]) . In the following, overdot denotes a r time derivative, and y A = y — z A . 

Noticing that the body zone radii are constant , i.e., R A = 0, we have 

P'a = P>a + Qa + ^, (18) 
4 = \(*ft + e*^pj + vM + le- 2 g-, (19) 

+ ^ + e 2 f» + ^, (20) 
Zf = -Af - A { } j)k , (21) 



where 



Mi ee 2e 4 / d*a A a\4; zi (22) 



gW = e -4 / rf5m (AJJ - w JAy z ) y*^, (23) 
R K / 3 = 6" 4 / rf5 m (A& - tff A&) yjy^, (24) 



and 



j T fc (*i) 

A*f* ee e 2 Jf «J + e\\jf + ifiM + e 4 ^— . (25) 

The symbol [ ] (or ( )) attached to some of the indexes denotes anti-symmetrization (or sym- 
metrization) on the indexes between the brackets. M A is the spin of the star A which is set to 
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be zero for the purpose of the current paper. Eq. ({TBI) gives a momentum-velocity relation. 
Thus our momentum-velocity relation is a direct analogue of the Newtonian momentum- 
velocity relation 60j. In general, we have 

rK,i 



j(K,i) + 21 j(Ki-x[ki)([ 



l + l 

z {mi + 



'A 



21 

7+T 



r(Kl-ilh)i}j 



+ z 



(Kij)i 



T+~L 



and 



J 



(Kii) 
A 



,dl< 



l + l dr 



-2lQ Kl 



y(Kii)j , 7 



(Kij)i 
A 



,2 (« rA"i)j 



Z + l 



A J 

2 



4 dJ A 



Kdij) 



l + l 



-21+2 -ftKiiij) 



(26) 
(27) 

(28) 
(29) 



/ + 1 dr 

Because the energy of the star, PJ, is not the mass of the star, we shall find the relation 
between these two in Sec. I VI I We may set the dipole moment D\ to be zero to define 
the center of the mass of the star z A , but we will later choose a different value for D l A for 



convenience. As for the multipole moments I A \ j^ 1 - 1 ^ 11 ^ anc i ^-'-m-wi^ we ma y equate 
those to zero, as it is our aim to derive equations of motion for a spherically symmetric star. 
However, there is a subtlety here, which we will explain in the next section. 



r (.?Q-i[fcl)i] 



AKi-i[ki)i]j 



2. A Spherically Symmetric Star and the Star's Multipole Moments 

As mentioned in the introduction, this paper studies a binary consisting of two spherically 
symmetric compact stars. In other words, all the multipole moments of the star defined in 
an appropriate reference coordinate where effects of its orbital motion and the companion 
star are removed (modulo, namely, the tidal effect) vanish. We adopt the generalized Fermi 



normal coordinate (GFC) 61] as the reference coordinate. 

We have defined the multipole moments of the star A as Eqs. f fT5l) . (TIB"]) , and (IT7|) . Those 
are defined on the r = constant three-surface and differ from ones defined in the GFC. Then 
a question specific to our formalism is if the differences between the multipole moments 
defined in Eqs. (fT5|) . ([TBI) , and (IT7|) and the those in the GFC give purely monopole terms. 
Roughly speaking, a body that is spherical in its GFC may acquire apparent multipole 
moments due to, for example, the Lorentz contraction. 

At the 3 PN order, this problem was addressed in the appendix C of Paper III and the 
differences were mainly attributed to the shape of the body zone. The body zone B A which 

16 



is spherical in the near zone coordinate (NZC) is not spherical in the GFC because of a 
kinematic effect (Lorentz contraction). In fact, the difference between the GFC quadrupole 
moment and the NZC one was found to contain monopole terms. 

At the 3.5 PN order, no NZC multipole moment was found to contain purely monopole 
terms. We shall be back to this point later in Sec. IV A 21 We equate all those NZC multipoles 
if (I > 3), jf l -^ [kl)t] (/ > 1) and zf l ~ l[mi (I > 1) to zero from now on in this paper. 



3. NZ/B Contribution 

About the NZ/B contribution, since the integrand A^ z = —gt^ L + x^ vafi ' , a p is at least 
quadratic in the small deviation field h^ u , we make the post-Newtonian expansion in the 
integrand. Then, basically, with the help of a (super-)potential g(x) which satisfies Ag(x) = 
f(x), A = did 1 denoting Laplacian, we have for each integral (, e.g., n = term in Eq. ( fTTT) ) 



d 3 y I® = -4ng(x) + I dS k 

NZ/B \ x V\ Jd(NZ/B) 



1 9g{y) d f 1 



x — y | dy k dy k \ \x — y 



(30) 



Eq. (!30~!) can be proved without using a Dirac delta functional (see Appendix B of Paper 
III). For n > 1 terms in Eq. fTTTl) . we use appropriate (super-)potentials many times to 
convert all the volume integrals into surface integrals and bulk terms ("— 47tg(x)") [67| . 

Finding the super-potentials is one of the most formidable task especially when we proceed 
to a high post-Newtonian order. Fortunately, at the 3.5 PN order, all the required super- 



potentials are available (See [34|, |38| and Sec. IV Bl below). 



C. General Form of the Equations of Motion 

From the definition of the four momentum (Eq. ([15"]) with / = 0) and the conservation 
law Eq. ([6]) we obtain an evolution equation for the four momentum; 

^ = -6" 4 / rf^A^ + 6- 4 ^/ dS k K^ z . (31) 

" r JdB A JdB A 

We have defined the four momentum by a volume integral of K r ^ z . Because x TUal3 ,ai3 = 
X TI/ai ,ai, we can transform the volume integral of X TUal3 ,a/3 info a surface integral form and 
can evaluate it explicitly in terms of the mass, velocity and the orbital separation. As a 
result, it is straightforward to see that x P ar t °f Eq. (!3~T!) is a trivial identity. The identity 
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to the 3 PN order was shown in Appendix E of Paper III by an explicit calculation. This 
observation implies that equations of motion can be derived from part of Eq. fl3Tj) . We 
thus define the part of the four momentum, dipole moment, and Q\ integral as 

P£ e = e 2 [ d 3 a A Q% z , (32) 
Jb a 

2 / J3, 



D\ e = e 2 / d 6 a A at A e T J z , (33) 
Jb a 

Qab = W dS k (6# z - «ley z ) y\. (34) 

J OB A 

Correspondingly, we split the momentum-velocity relation Eq. ffl8l) and the evolution equa- 
tion for the four momentum Eq. (I3TI) into the part and the x part. Note that the x P ar t 
of the multipole moments and the Q A ' 1 and R A 1 ^ integrals still in principle affect the equa- 
tions of motion through the field, as ,a/3 emerges as the difference between the curved 
spacetime d'Alembertian and the flat spacetime one and thus should affect the gravitational 
field. Indeed the x P ar t of the star's energy shown in Appendix |A] affects the acceleration 
at the 3.5 PN order. 
The part of (HI) is 



^ = -e" 4 / dS k Q% + e~\ k A l dS k e T J z . (35) 
dT JdB A Job A 

Substituting the part of the momentum-velocity relation into the spatial components of 
Eq. (1351) . we obtain the general form of the equations of motion for the star A; 

= -e~ 4 / dS k Q% z + e~%\ I dS k 0™ z 



dr 



dB A JOB A 



+e~ 4 v\ ( I dS k Q k N T z -v\<f> dS k Q T N T Z ) 

\JdB A JBBa J 



dr dr 2 

All the right hand side terms in Eq. ( |36l) except for the dipole moment are expressed as 



surface integrals. We can specify the value of D l AQ freely to determine the representative 
point z a (t) of the star A. 

In Eq. ( !36l) . P Ae rather than the mass of the star A appears. Hence we have to derive 
a relation between the mass and P A q- We shall derive the relation by solving the temporal 
component of the evolution equation (1351) functionally. In fact, at the lowest order, we have 
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shown in Paper II that 



dPl 
dr 




(37) 



Then we define the mass of the star A as the integrating constant of this equation; 



ttla is the ADM mass that the star A had if the star A were isolated. We took e zero 
limit in Eq. (1381) to ensure that the mass defined above does not include the effect of the 
companion star and the orbital motion of the star itself. Some subtleties about this definition 
were discussed in Paper II. By definition itla is constant. The procedure that we solve the 
evolution equation of P t A q and obtain the mass energy relation is achieved up to the 3.5 PN 
order successfully and the result will be shown in Sec. PVTl 

D. On the Arbitrary Constant Ra 

Our final remark in this section is on the two arbitrary constants Ra- Since we introduce 
the body zones by hand, the arbitrary body zone radii Ra seem to appear in the metric, the 
multipole moments of the stars, and the equations of motion. Paper II has proved that the 
surface integrals in the general equations of motion Eq. (I36I) do not depend on Ra through 
any order of the post-Newtonian iteration. The appendix D of Paper III explained that the 
field and the multipole moments are independent of €Ra- 

Practically, those two accounts justify that we safely discard all the cRa dependent terms 
except for logarithms of €Ra that appear in the course of computations. We keep lnei?^ 
dependent terms to make the arguments of the logarithms adimensional. 

We here emphasize that we discard the cRa dependent terms in the field first and then 
evaluate the surface integrals in the general form of the equations of motion using the field 
independent of €Ra- We then discard the sRa dependent terms arising in the computation 
of the surface integrals. The details of this procedure were explained in Paper III. 

IV. STRUCTURE OF THE 3.5 PN EQUATIONS OF MOTION 

In the following sections, we shall derive an acceleration for two spherical compact stars 
through the third and a half post-Newtonian accuracy. For this purpose, we evaluate the 



m A = \imPX . 



(38) 



19 



surface integrals in Eqs. (|35|) and (|36|) to the appropriate order. As the mass of the star is 
0(e 2 ), the Newtonian force appears as e 2 correction to the lowest order equations of motion 
(m A dv l A /dT = 0). The 3.5 PN order correction, or (f/c) 7 correction to the Newtonian force 
appears at 0(e 7 ) and hence the equations that we have to evaluate to derive an evolution 
equation for the energy and the equations of motion are 



dP\e \ _ ( dPie \ , £ 7 



d T ' <3.5PN V dr /<:;!'. N 

. dv\ \ ( dv\ , 

m x [ — J = m 1 I — I + e 

CtT / <3.5PN V " T / <3PAf 



* dS k u Q t nz + Vi f dS k u Q t nz 

JdBx JdB x 



(39) 



ds k u^Nz+n 4> ds k u e^ 

aBi 



V aT J 3.5PN V Ur / 3.5PN 

rfr dr 2 ' 1 J 

where for an equation or a quantity /, (f)< n pN and (f) n PN denote / up to the n PN order 
inclusively and / at the n PN order, respectively. < n f and n f on the other hand denote 
an equation or a quantity / up to 0(e n ) and at 0(e n ), respectively. In Paper II, we found 
Q Ae = 0(e 6 ). It should be understood that in the second line of Eq. (HDjl . the acceleration 
dv\/dr should be replaced by the acceleration of an appropriate order lower than the 3.5 
PN order. Henceforth, we call Eq. (140]) the general form of the 3.5 PN equations of motion. 
The explicit forms of the integrands = n[— gfc^A ( on ®Ba) are 

16tt u ey z = - 7 -,h T \ k 9 h^' k + • • • , (41) 
IQtiu&nz = ^h T \ k 9 h T[k ' i] + --- , (42) 
16tt n Q% z = l -(5\5h + 5* A - 8^5 M ) { 4 h rT '\ u h TT ' 1 + 9 h m J + 4 9 h r \ T ) + 8 4 h T m > k 9 h T ^} 
+ 2 4 h T \ k 9 h T ^ + 2 A h T \ k 9 h T ^ + ■■■ , (43) 

See the appendix [B] for the complete expressions. The field components up to the 2.5 PN 
order inclusively, < 9 h TT , <7h T \ and , are listed in Paper II. In addition to those, to 
derive the 3.5 PN mass-energy relation and the 3.5 PN momentum-velocity relation, we 
have to derive 9 h Tl . To derive the 3 PN equations of motion, we further need \\h TT + 9 h k k . 
We do not need to use the 3 PN field 8 h T \ 8 h^, and ioh TT for the purpose of this paper. 
Up to the 2.5 PN order inclusively, the super-potentials required to compute the field could 
2s| . At the 3 PN order, it is quite difficult to complete the required super- 



be found 
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potentials. We took another method to overcome this problem in Paper III. Fortunately, 
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we can find the super-potentials at the 3.5 PN order and hence derive the 3.5 PN field in a 
closed form. The next section shows our method of the derivation of the 3.5 PN field. 

V. 3.5 PN GRAVITATIONAL FIELD IN HARMONIC COORDINATES 

As written in the previous section, our derivation of the 3.5 PN equations of motion 
requires gh Tl and \\h TT + gh k k- We start with the body zone contribution. 

A. 3.5 PN Body Zone Contribution for a Spherically Symmetric Star 

This paper studies equations of motion for a spherically symmetric star and hence discard 
all the GFC multipole moments of the stars. As stated in Sec. IHIB 21 and will be shown in 
Sec. IVA21 up to the 3.5 PN order we can safely discard the NZC multipole moments I A \ 
jiKt-ilhW and giK^AhW except for the NZC quadrupole moments PJ. Using Eqs. JTSJ) to 

(|29|) . the body zone contribution for a spherically symmetric star becomes 





(44) 




5 




(45) 



P 



<r„,i „J 



5 



Q 




+ 5> 



n=4 



5 



.ii 




,klm(ij) <klm> 
'A 'A 



+ 



n-n-A 



, r A l r>kji , r>kij _ Tjijk\ 
~T 3 {n-n-A n-H-A n-^A I 



A 




(46) 
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In the above equations, D\ may be specified freely to define the star's representative point 
z\. Although we explicitly wrote Q A ' 1 and R A 1 ^ integrals above, the next section shall show 
those do not affect the 3.5 PN acceleration. 



1. Q K / and R K / j 

The surface integrals Q A li and did contribute to the field and the equations of 

motion at the 3 PN order (Paper III). 

At the 3.5 PN order, we compute Q A 11 and up to 0(e 5 ) to derive the 3.5 PN 

field. Then the integrands of the surface integrals in Eqs. (!23l) and fl24l) are gA^ z . Now 
suppose those integrands behaves near the star A as l/r A with p some positive integer. Then 
Qa' 1 = 0((€Ra) 1+3 ^ p ) and R A ' % ^ = 0((€Ra) 1+4 ~ p )- Because we neglect any explicit terms 
that depend on (eiiU) 9 (q- non-zero integer. See Sec. IIIIDI and Sec. Ill E of Paper III), 
we compute only Q A l% for < I < p — 3 and R A ll: * for < I < p — 4. Inspection shows 
that p = 6 for gA^f z and p = 7 for gA^ z . For those degree I of the collective indexes K h by 
explicitly computing those integrals in Eqs. (1231) and fl24|) . we found that Q A 11 and R A [l:j do 
not contribute to the field nor equations of motion at the 3.5 PN order. 



2. A Spherical Body and Multipole Moments 

As promised, this section explains the reason we can safely discard the NZC multipole 
moments except for the NZC quadrupole moments through the 3.5 PN order. Denoting 
the NZC (GZC) quadrupole moment as I ANZC (-Oigfc)' (^14) of Paper III gave an 
expression of the difference between these two moments; 

= 'Xnzc - 'Zgfc = e~ 8 (B\t)v a - A«(r)) / dS m y A y A y%A T N T z , (47) 

JdB A 

for a star which is spherically symmetric in the GFC. The coefficients A 1; >(t) and B l (r) can 



be read off from 61| (see also Eq. (C16) and (C17) of Paper III) and are 



B\t)v{ = eVi + 0(e 4 ), (48) 
AV{t) = e 2 (\v[v{ - ^5 A + 0(e 4 ), (49) 



.2 1 1 r 12 

for the star A = 1 and fi2 = Z\ — zi = r<i — t\ . Note that there is no e 3 term in A^{r) nor 
B 1 (t) because there is no 0.5 PN field in the harmonic gauge (i.e., ^h TT = 0) as was shown 
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in the appendix E of Paper II. The surface integral in Eq. (1471) can be expanded in e as 
e- 8 <L dS m y A y A y^A T ^ z = f dS m y A y j A y^A^ z + e f dS m y A y j A y^ 9 A T ^ z 

JdB A JdB A JdB A 

+ 0(e 2 ). (50) 

The first integral in the right hand side of the equality is found to be non-zero and contribute 
to the 3 PN metric ioh TT and hence the 3 PN equations of motion. The second integral 
could in principle affect the 3.5 PN field, as the coefficient e in front of it indicates. So let us 
evaluate the second integral. Its integrand is gAy z = — iJi TT ^h % i which behaves near the 
star A = 1 as gAy z ~ 3r\r{/rl — 5^ /rf. The surface integral then results in 

/ dS m y\y\y^A^ z = 0(e 2 R 2 A ), (51) 

JdB A 

and thus there is no monopole term in SI A at the 3.5 PN order. Similarly, using Eq. (C14) 
of Paper III, we can check that none of the NZC multipole moments hide any monopole 
terms at the 3.5 PN order for a star spherically symmetric in the GFC 



B. 3.5 PN N/B Field 



This section explains how we derive the 3.5 PN N/B field components 9h T ^ z ^ B and 
h^-jvz/b + Q^NZ/Bk- The integrals for the latter which we evaluate may be written as 

nh T N Z/B (T, x) + 9h k NZ/Bk (r, x) 
= 4 V ( ~ 1)n dn ! ^3 ii-nA]y r z(r,y) + 9-A fc (^y) 

- 4 / dS j l0 A T J z (r,y) - [ d 3 y 8 A k NZk (r,y). (52) 

JdNZ ° T JNZ/B 

Because only a spatial derivative of yJ^zib ~^^NZ/B k a PP ears i n n[ — 9^ll\ as shown in Eq. 
f )43|) . the second and the third term in the right hand of the equality in the above equation 
does not contribute to the 3.5 PN equations of motion. Similarly, 9h^ z , B may be written as 



9^z/b(t, x ) = 4 / d V A ^ z{T f - 4 / dS j s A% z (t, y) 

JNZ/B F ~~ V\ JdNZ 

f d 3 y 6 A]$ z (r,y)\x-y\ 3 . (53) 

JNZ/B 



'NZ/B 

2 d 3 



3<9r 3 
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Our task now is to find super-potentials / that satisfies, e.g., nA^ + gA^fc = Af. With the 
super-potentials, we use Eq. f l30|) and find explicit expressions of the 3.5 PN field components 
in closed forms. For the higher order retarded expansion terms, for example for the second 
retarded expansion term of Eq. (|52j) (the n = 2 term of the first term in Eq. fl5"2"j) ) , we use 
super-super-potential f(y) satisfying gA^ z + jA k NZk = A 2 f. Then the volume integral may 
be converted into the bulk term and the surface integral terms as 

/ d 3 y\x-y\A 2 f(y) 
Jnz/b 

= Snf(x) 

yk _ x k 2 2(y k — x k 



+ f dS k 

'd(NZ/B) 



x - y\d k Af(y) - !L__^Af(y) + ^^d k f(y) + » ' f{y) 
\x — y\ \x — y\ \x — yy 



The source terms for which we need to find particular solutions of Poisson equations have 
the following spatial coordinate dependence. 
111111 

— fi ' K ' ~~ 4 ' 3 ' 2 ' ' 1 ' 

ry 1 - 1 ry'-J ry^ y <-* ry '** ry 
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(54) 



and 1 «-> 2, i.e., like l/r| for (Specific members of the list depends on how one simplifies 
the expressions of the source terms of the Einstein equations. For instance, one can always 
erase r\ by using r l 2 = r\ 2 + r i-) 

Our method to derive the super-potentials are heuristic; there are few guidelines available 
to find the required super-potentials. We proceed as follows. First, we convert all the 
tensorial sources into scalars with spatial derivatives. For example, nA T ^ z + $A k NZ k include 
the following term 

Z 2 \'2 



. (n • ^)(f 2 • V)(r 12 ■ V) = -2A mi m 2 2v\V^r 12 ■ V) — — ^ . (55) 



(Here and henceforth, it should be understood that in general "scalars" can have tensorial 
indexes carried by va and f 12 , but do not have those by ta-) 
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Second, we find the particular solutions for Poisson equations with the scalars as sources 
using a formula A(f(x)g(x)) = g(x)Af(x) + 2V/(x) ■ Vg(x) + f(x)Ag(x) valid in NZ/B. 
We also use super-potential chains such as; 

(-3,-2) Aij^ Rf (-5,-2) 



where 



| A 22 I A22 

2j(-3,-4) ^ 12 j(-5-4) 

/( -3,-2) 



ln[^ 

n 



an d /K") satisfies A j(m,n) = r m r n = / Q^J Q Za , 

easy to find a particular solution, 



For example, for Eq. (155]) . it is 



lit 
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1 d d 



if 



where 
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(4,-3) 



2^,2 



rfr 



4r 9 



1' 12 
2ro 



12 <9^ dz, 



+ r^r 2 + 2r^ 2 r 2 



(4,-3) 



1r\ rf 2 

mr 2 . 

3 r 2 



(56) 



Following the method described above, we could find all the required particular solu- 
tions and using algebraic computation codes written in Mathematica we have derived the 
expressions of gh Tl and \\h TT + gh k k in closed forms. However we do not reproduce those in 
this paper because the numbers of terms in the expressions are huge (~ 800 and ~ 1000, 
respectively. These numbers depend on a particular simplification one makes on those ex- 
pressions) . 

Incidentally, we have checked (a part of) the 3.5 PN harmonic condition <g/i" M ,^ = 0. 
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VI. 3.5 PN MASS-ENERGY RELATION IN HARMONIC COORDINATES 



By evaluating the surface integrals in the evolution equation of the star's energy Eq. (1391) . 
one may obtain the time derivative of the mass of the star 1 as 
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5 



Paper III gives the explicit expression of (dP^ e /dr) <3PN . 
We can integrate Eq. (ISTj) functionally as 

7 

Pie = m 1 ^e\r 1 + 0(e 8 ). 

fc=0 



(57) 



(58) 



The a up to 3 PN order are given in Paper III. The new result we show in this paper is 
7!" 'a; 



7I1 



^m 1 2 m 2 (ni2 • V) ^m\m 2 {ji\ 2 ■ V)V 2 I6mim 2 2 (ni 2 ■ V) 



15ri- 3 



12 



15ri- 2 



12 



(59) 



where V = V\ — v 2 . The mass-energy relation for the x P ar t up to the 3.5 PN order is given 
in the appendix [A] Eqs. ( 1581) and ( 1A2I) give the 3.5 PN order mass-energy relation in our 
formalism. 

From the definition of P A q, we expect it equal to y/—gm A u T A where u r A is the time compo- 



nent of the 4-velocity of the star A normalized as g^ v u A u A = — e with u\ = u T A v 



~9U A 
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at 3.5 PN order can be written in terms of the deviation field as 



-9U A 



<6l 



+ e' 



4 ! 



,h TT - 



7 h\ - 7 h n v\ + - 7 h TT v 2 A 



+ ^hijV A v A - - 5 h k k v A + ^ b h k ki h 7 



16 

0{e% 



(60) 



This expression should be somehow evaluated at the star A. However, because the deviation 
field h^ u diverges at the star in the point particle description, the above expectation is not 
trivial. As in Paper I, II, and III, we checked that up to the 3.5 PN order inclusively a 
relation 



P A0 =m A [ 



-9u\rz\ 



(61) 



holds where [f] A xt means that we regularize the quantity f at the star A by the Hadamard's 
Partie Finie (see e.g. 46j) or whatever regularization which gives the same result. We here 
emphasize that we have never assumed this "natural" relation in advance. The relation Eq. 
flBTj) has been derived by solving the evolution equation for P Ae functionally. 

For the 3.5 PN mass-energy relation Eqs. fl58l) . one can check that only the field com- 
ponents up to the 2.5 PN order <gh TT and < 7 h tu appear in the expression of y/—gu T A in 
Eq. fIBTl) above. We apply a regularization on those components in the right hand side of 
Eq. ( JBTT) . The "naturality" of Eq. ( JBTT) thus supports a use of the Hadamard Partie Finie 
regularization in the literature up to the 2.5 PN order 28]. 



VII. 3.5 PN MOMENTUM- VELOCITY RELATION IN HARMONIC COORDI- 
NATES 

From Eq. ( fl8l) . an explicit expression of the momentum- velocity relation is obtained by 
evaluating the Q AQ integral up to 0(e 7 ) 

Qab = <sQab + e 7 / dS k (n[-9t T L k L } - n[-9t T L T L }v k A ) y\. (62) 

JdB A 

The computation is straightforward and we found 7 Q % A ® = 0. From the 3 PN accurate Q l Ae 
calculation, we have a momentum-velocity relation. 

Pie = Phv\ ~ e 6 ± Qmfai) + e 2 ^ + 0( e 8 ), (63) 
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where a\ is the acceleration of the star 1 and should be replaced by the 3 PN order expression 
of the acceleration. As in Paper III, we define the representative points of the stars z\ by 
choosing 



D\q{t) = e^m\a l A - e 4 ^-m A (j A \n (-^~ 



(64) 



The first term in Eq. (j64p makes the three momentum proportional to v A at the 3 PN order. 
The second term gauges away the logarithmic terms from the 3 PN equations of motion 36]. 
In any case, our choice of D Ae affects only the 3 PN order correction and P{ B = P[qv{ holds 
as long as we are concerned with the 3.5 PN order correction to the equations of motion. 

As the second term in Eq. (|64|) depends on the arbitrary parameter eR A , one may suspect 
that our equations of motion lose their predictive power on the binary dynamics. This is not 
the case. When we define the star's representative point by D Ae = instead of Eq. (JMl) . 



logarithmic terms that depend on eR A appear in equations of motion [36(. But those terms 
are mere gauge. Indeed, it can be shown that an observable such as the orbital energy does 
not depend on eRA when that observable is written in terms of gauge independent variables 
(such as the gravitational wave frequency) 29j, |35 |. 



The shift in the world line induced by the second term in Eq 



is equivalent to the 

gauge transformation used in the works by Blanchet and Faye |29j, l32|. However, the fact 
that we adopt Eq. (1641) does not mean that we use a regularization in any sense. The 
generalized Hadamard Partie Finie or the dimensional regularization is nothing to do with 
any mere shift of the world line. In fact, the undetermined coefficient A associated with 
their use of the generalized Hadamard Partie Finie cannot be gauged away by any shift of 



the world line 



29 



22]. 
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VIII. 3.5 PN EQUATIONS OF MOTION IN HARMONIC COORDINATES 



With the 3.5 PN field at hand, we evaluate the surface integrals in the 3.5 PN general 
equations of motion Eq. (T4"0l . A tedious but straightforward calculation results in 



(mial)< 3 .6PN = (miai)<3PN 
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where the acceleration up to the 3 PN order is given in Paper III. 

Eq. ( 1651) is in perfect agreement with the previous works in harmonic coordinates 



(65) 



the result in the ADMTT coordinate 



10 
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results from the energy balance argument 



38|, 



15 ] by a suitable gauge transformation, and also the 



16,HQ. 



We have used the local conservation 



law of the stress energy tensor of the matter and the gravitational field and the surface 
integral approach to derive our 3.5 PN equations of motion. We have not a priori assumed 
that the star follows a geodesic in any sense. The strong field point particle limit enables us 
to realize a point particle with strong internal gravity without using a Dirac delta functional. 
Nissanke et al. [381] assumed that a star follows a geodesic regularized by the Hadamard 
Parti Finie regularization (or any other regularization method that gives the same result, 
such as the dimensional regularization). Thereby, the perfect agreement between our present 
work and that work |38] confirms that a self-gravitating star follows the regularized geodesic 
at least up to the 3.5 PN order inclusively. 
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APPENDIX A: X PART 



This section shows the functional expressions of P\ x in terms of tua, v a , V 1 = v{ — v 2 , 



and r\ 2 . Here we defined as 



^a x - e 



1 ' rfW^V 

'B A 



(Al) 



By the definition of x^ val3 



a/3, 



167i X TTa(: \ a p = {h Tk h Tl - h TT h kl ) 



kl; 



and thus we can obtain the functional expressions of P£ by evaluating surface integrals 
using the Gauss's law. In fact, up to the 3.5 PN order, the definition of Eq. flAll) gives 



pr -Apr, • V) 6 pT 
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45 9 5 

A6v 2 2 (f 12 ■ v 2 ) 196(ni2 • v\){v\ ■ v 2 ) 80(n 12 • v 2 )(vi ■ v 2 ) 

45 + 9 

20(ni2 • v{){n l2 ■ v 2 ) 2 + 8(n a2 • v 2 f 



45 



9 

78(ni2 ■ vi) 2 (n 12 ■ v 2 ) 



(A2) 



The explicit expressions for the 2 PN and 3 PN terms, 4-Pi T x , and &P{ x i are given in Paper 
III. The 3.5 PN \ part of the star's energy iP\ x does contribute the 3.5 PN field \\h TT (See 
Eq. ( 1441) ) and affects the 3.5 PN equations of motion. 
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APPENDIX B: 3.5 PN LANDAU-LIFSHITZ PSEUDO-TENSOR 



This section lists the components of the Landau-Lifshitz Pseudo-tensor at 0(e u ) in our 
ordering which are necessary to compute the 3.5 PN equations of motion in our formalism. 
See the Paper II for <g[— lGngt^] and the Paper III for 10 [— IFmgt^j]. Note that all the 
divergence such as h^ k ^ in the paper II should be replaced by — /i MT jT in consistency with 
the following results. This is simply because it is practically much easier to evaluate h tJ,T jT 
than —h^ k tk . 

The Landau-Lifshitz pseudo-tensor 65j in terms of which satisfies the harmonic 
condition is 

(-16^)^ = g a ^ 5 h^, 7 h^, s + \g^g a ph a \ s h^ n - 2g aP giW jS h s ? 

+ \ (g m g u0 - (g^g* ~ ^gsc) h^, a h s <^. (Bl) 

We expand the deviation field h^ u in a power series of e; 

n=0 

Paper II showed the lowest order of the field is e 4 . Using this expansion, we expand t^ L in 
e. The point to note is that 1) we raise or lower indexes with the flat metric 1]^ and rj^, 2) 
■rfT = - e 2 and r] TT = -e~ 2 , 3) b h l \ k = and 7 h T \ k = 0, 4) 5 /i T ^ = 0. 

In the following, repeated alphabetical indexes (i.e. excluding r) must be summed. In- 
dexes in round (square) brackets, (...) (or [...]) ,means (anti)-symmetrization on the indexes, 
and indexes between vertical bars, |...|, are excluded from (anti)-symmetrization. 

1,3 , 7 , 

ii [ — S^ll] = 7 ih TT tT $h \,T — 7 4^ rT ,r 7^- tt ,t — jh T >T 4,h TT >k + ^ jh TT Ji TT ' /Ji TT \ k 
4 4 o 

7 7 

+ -^hki Ah TT ' k 4h TT ' 1 — — 4^ rr ,fc 9h TT,k + §h k i T 4h Tk ' L + 2 4 /i T k , 7 h r( - k '^ 
o 4 

+ \ 4h TT ' k 7 h l l>k . (B2) 
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+ I 7 h kl 5 ij 4 h TT ' k 4 h TT ' 1 

+ I rh ij 4 K\ k 4 h^ k - \ 7 h k ^ 4 h^ 4 h^ k 
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